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Abstract 


In social and behavioral sciences, data are typically not normally distributed, which can 
invalidate hypothesis testing and lead to unreliable results when being analyzed by 
methods developed for normal data. The existing methods of generating multivariate 
non-normal data typically create data according to specific univariate marginal 
measures such as the univariate skewness and kurtosis, but not multivariate measures 
such as Mardia’s skewness and kurtosis. In this study, we propose a new method of 
generating multivariate non-normal data with given multivariate skewness and kurtosis. 
Our approach allows researchers to better control their simulation designs in evaluating 
the influence of multivariate non-normality. 

Keywords: multivariate non-normal data, multivariate skewness, multivariate 


kurtosis, random number generation 
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A Method of Generating Multivariate Non-normal Random Numbers with Desired 


Multivariate Skewness and Kurtosis 
Introduction 


In social and behavioral sciences, the normality of data is assumed in most 
statistical methods. Nonetheless, data are rarely normally distributed in practice. 
Therefore, the statistical inferences may not be valid, and the results may not be 
reliable any more when procedures developed for normal data are used to analyze 
non-normal data (Cain, Zhang, & Yuan, 2017; Micceri, 1989). Many studies in the 
literature investigated the consequences of the violation of the normality assumption 
and proposed some alternative procedures to analyze non-normal data. For instance, 
Bradley (1980) showed that robustness of statistical procedures such as the classical Z, 
t, and F' tests suffered from the non-normality of data. Non-parametric tests and 
procedures have won appreciation of researchers because they do not rely on the data 
distribution and therefore, the violation of normality does not directly disqualify data 
analysis (Hollander & Wolfe, 2015). 

In the literature, discussions on non-normality mainly focus on the univariate 
case; whereas the consequences of deviation from the multivariate normality is less 
explored. However, the analysis of multivariate data is routinely conducted in social 
and behavioral sciences research. Therefore, it is important to understand the influence 
of the multivariate non-normality on the multivariate analysis, which can be done 
through Monte Carlo simulations. To conduct such simulations, one needs to generate 
multivariate data with the control of the degree of non-normality. In the literature, 
most non-normal data generators are developed for univariate data, such as the 
third-order polynomial power method (the power method; Fleishman, 1978), the 
fifth-order polynomial method (Headrick, 2002), and the g-h distribution method (Field 
& Genton, 2012). 

The existing methods typically generate multivariate data according to specific 
univariate marginal measures such as the univariate skewness and kurtosis, but not 


multivariate measures such as Mardia’s (1970) skewness and kurtosis. For example, the 
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widely-used simulation method proposed by Vale and Maurelli (VM; 1983) was built on 
Fleishman’s (1978) polynomial approach. In addition to generating data for each 
variable with specific first four moments, their method also controls for a correlation 
matrix that allows researchers to have a desired multivariate data structure. This 
method is very popular in the moment-based modeling area, such as structural equation 
modeling (SEM). However, some researchers have questioned the generalization of this 
method. Foldnes and Grgnneberg (2015) derived the mathematical distribution of the 
VM approach and showed that even though the approach could generate multivariate 
data with user-specified marginal skewness and kurtosis, the generated data might not 
be truly multivariate non-normal. Astivia and Zumbo (2015) have shown that the Vale 
and Maurelli method has downward bias. In one of their later papers, Astivia and 
Zumbo (2018) also found the multiplicity solution issue of the Fleishman’s polynomial 
related method, which means that there are multiple possible solutions for the 
polynomial coefficients (a, b, c, and d). This issue might lead to the difference in the 
analysis even with the same inputs. To remedy the drawback, researchers have 
developed other methods. Mair, Satorra, and Bentler (2012), for example, introduced a 
multivariate approach based on copulas that could also generate data with a 
pre-specified variance-covariance matrix. Foldnes and Olsson (2016) presented a method 
using linear combinations of independent generator variables. Additionally, Lee and 
Kaplan (2018) developed a generator for the multivariate ordinal data based on entropy 


procedures. 


Despite their usefulness, none of these methods allows the direct control of the 
multivariate non-normality measures. Multivariate skewness and kurtosis have been 
shown to directly impact statistical analysis. For example, Yuan, Bentler, and Zhang 
(2005) noted that a robust procedure might be necessary for reliable SEM inferences 
when a sample has a large multivariate kurtosis. More recently, Cain et al. (2017) 
conducted a meta-analysis study on the multivariate non-normality of the data used in 
254 published studies and found that the Type I error rates of testing the model fit 


were remarkably higher in factor analysis when the multivariate normality was violated. 


GENERATING MULTIVARIATE NON-NORMAL DATA 5 


Generating multivariate non-normal data with desired multivariate measures is the first 
step to understand the type and severity of non-normality. This is because it relates to 


both multivariate skewness and kurtosis on analysis procedures. 


Generating multivariate non-normal random data requires the understanding of 
the definition of the non-normality and the relationship between univariate and 
multivariate data. Mardia (1970) introduced the measures of population multivariate 
skewness and kurtosis as the natural extension of the univariate ones. In the univariate 
case, with non-zero skewness, the distribution is asymmetry. When the excess kurtosis 
is not 0 (excess kurtosis equals to kurtosis minus 3), the distribution density function is 
different from a normal distribution. Similarly, Mardia’s multivariate kurtosis indicates 
whether the tails are heavy or light in comparison to those of the multivariate normal 
distribution (DeCarlo, 1997). On the other hand, the Mardia’s skewness is still a 
measure of symmetry, but cannot take negative values. Higher values indicate severer 


asymmetry. 


To date, there is no available method for researchers to directly specify both 
multivariate skewness and kurtosis for multivariate non-normal data generation. To fill 
the gap, we introduce a new method of generating multivariate non-normal data with 
specific multivariate measures. This approach allows researchers to better control their 
simulation design in evaluating the influence of the multivariate non-normality. More 
over, this technique will allow for a better understanding of the relationship between 


multivariate non-normality and the marginal univariate non-normality. 


The rest of the paper is organized as follows. We first propose a new generating 
method and introduce an R package for the implementation of the method. We then 
present a simulation study and the results with various conditions. We conclude the 


study with a summary of our method. 
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Method 
Data Model 


To generate the non-normal data, we specify the following data model. We use a 


vector x of p variables as, 


x=TrA€, (1) 
and each marginal x; as, 
qd 
m= 1 Sails, (2) 
j= 


where € = (&1,...,€,) is a vector containing g independent random variables. Each of the 
variable €; has first four ordered moments E(€;) = 0, E(€7) = 1, E(&?), and E(€5). 

A = (a;;) is a p X q matrix of rank p (p < q), and AA’ = © = cov(x). And risa 
random variable, which is independent of €, with the first four ordered moments 

E(r), E(r?), E(r?), and E(r*). 

The ordered moments are a set of quantitative measures describing the shape of a 
distribution. When the ordered moments are normalized, they become the standardized 
moments (or central moments). Skewness and kurtosis are the third and fourth 
standardized moments. The ordered and standardized moments can convert to each 
other as long as mean and variance are provided. 

According to the definition of Mardia (1970), the population multivariate 
skewness ( (3; ) and kurtosis ( (62 ) of x based on our model are computed as, 


B= B{{(x — w)'B'(y — WP} = Br) NEE (3) 


j= 


BR 


Ba = B{ l(c — wy E"x - ica -ye 


Using these formulas, population multivariate skewness and kurtosis can be calculated 


when univariate measures”, the ordered moments of r and €;, are given. The 


? In this paper, when it is in the multivariate setting, we refer the moments of x; as marginal measures. 
And the univariate measures are the moments of r and €;. The marginal skewness and kurtosis of x; are 
m(2i) = B(r®) OF a (G)/oif”, 

y2(a;) = E(r*)[ a aj, (E (€}) — 3)/o?, +3] (Yuan & Bentler, 1997). 


GENERATING MULTIVARIATE NON-NORMAL DATA ve 


standardized multivariate kurtosis formula, centering 32 by p(p +2), has been obtained 
in Yuan, Zhang, and Zhao (2017). However, the solution of univariate measures cannot 
be uniquely obtained based on these formulas from specified multivariate measures. 

Although the solution is not on an one-to-one basis, multiple solutions are 
available that sharing the same multivariate measures. Since we only care about the 
measures at the multivariate level other than the univariate level (or the marginal 
level), to remedy the lack of uniqueness of the solution, we establish one from 
multivariate to univariate by applying some constraints. 

First, we set r to be a constant 1 for convenience because it is only a scale factor. 
Second, the number of variables in € is set to be the same as the number of variables in 
X, so that p = q. Additionally, ; to € are set to be independent and identically 
distributed (i.i.d.). Thus, the 3rd and 4th ordered moments are the same for all §;, 
which are defined as E(€?) and E(€*). With these constraints, the multivariate 


skewness and kurtosis above become, 


Bt = E{{x— py ="'(y — w)P} = ABO YP, (5) 

By = E{[(x — w)'="(x — w))?} = pE(E*) + p(p — 1). (6) 
When multivariate measures (Gf and 33) and number of variables (p) are provided, 
E(€?) and E(€*) can be computed through Equations (5) and (6). Once we have the 
moments for all €; (i.e., E(€) = 0, E(€*) = 1, E(€°) and E(€*)), we can generate 
random numbers of € and then transform them to multivariate random numbers of x by 
Equation (1). 

We acknowledge that relaxing some of the constraints would give researchers more 
control of the univariate measures (e.g., allowing different univariate measures or setting 
different scaling factors). However, based on a small-scale simulation, we found they did 
not influence the behavior of the multivariate measures much (see Table 2 in the 
supplementary materials). Therefore, we use the constraints for convenience in this 
study. 

For each €;, we use a modified power method to generate random non-normal 


numbers. The widely-used power method is proposed by Fleishman (1978), which is to 
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generate non-normal data through a polynomial transformation 
Y=at+bZ4+cZ? +42’, (7) 


where Z comes from a standard normal distribution. With the information of the first 
four desired standardized moments (mean, variance, skewness 7;, and kurtosis 72) of Y, 
Fleishman (1978) derived four equations to obtain the coefficients a, b,c, and d through 
Newton’s method. 


In our study, Y is replaced by each €; as 
&=a+bZ+cZ? +dZ°. (8) 


Instead of including the standardized moments (skewness and kurtosis) as used in 
Fleishman’s method, we use the third and fourth ordered moments of €;. Therefore, the 


four equations of solving the coefficients are revised as 


a+c=0, (9) 

b* + 6bd + 2c? + 15d? —1=0, (10) 

72bcd + 6bc + 8c® + 270cd? — E(E°) = 0, (11) 

3b4 + 60b?c” + 60c* + 60b'd + 936bc7d + 630b7d? (12) 
+4500c?d? + 3780bd*® + 10395d* — E(€*) = 0. (13) 


With the value of a,b,c, and d, we first generate random numbers from the 
standard normal distribution to form the sample Z with size n. Then, the sample of €; 
is obtained by the polynomial transformation in Equation (8). Repeatedly sampling Z 
and conducting transformation for each €;, one gets the multivariate data with 
€ =(&1,...,€). The final step is to obtain x by applying the specific covariance matrix 
A to € following the data model in Equation (1) with r = 1. 


In summary, the following procedure can be used. 


1. With the user-specified multivariate skewness (7) and kurtosis (33) and the 
number of variables (p), calculate the third and fourth ordered moments of €; 


(f=1,2,-- ,p). 
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2. Generate the standardized €; by the modified power method to form €. 


3. Use the Cholesky decomposition to decompose the user-specified correlation 


matrix (or covariance matrix) to matrix A and multiply it to € (x = A&). 


Through this process, the generated data x will have the desired multivariate 
skewness and kurtosis. One shortcoming of applying the Cholesky decomposition 
approach is that, after the linear transformations, the population marginal measures of 
x (e.g., y, and 2) will be different from the original univariate measures of €. However, 
unlike the VM method, where the marginal measures are of interests, the focus of our 
method is the multivariate measures and the marginal measures are nuisance 
parameters. Therefore, our method does not require an intermediate correlation matrix 


and can apply the Cholesky decomposition directly. 


Limited Ranges of the Skewness and Kurtosis 


The power method cannot cover all the possible combinations of univariate 
skewness and kurtosis. This is because the method does not require the distribution of 
Y, and thus the moments cannot be analytically derived. However, the range 
relationship of univariate skewness (7,) and kurtosis (y2) of Y has been estimated 
through simulation (Luo, 2011). Based on that relationship, we derived the range 
relationship between the univariate €,’s third and fourth moments with our modified 
power method, which is 


E(€é*) > 1.641[E(€°)]? + 1.774. (14) 


Plugging it into the data model in Equation (1), the relationship of x;’s skewness and 


kurtosis is, 


1.641 
72> Scared ~ 1.226 at, +3, (15) 
ij ij 


which is restricted compared to the theoretical relationship of the general univariate 


skewness and kurtosis, 


m2 yt. (16) 
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Correspondingly, applying the inequality in Equation (13) to the multivariate skewness 
and kurtosis formulas in Equations (5) and (6), the relationship of multivariate 


skewness and kurtosis in our method can be derived as 


Bz > 1.64167 + p(p + 0.774). (17) 


R package 


An R package mnonr is developed based on our method to generate multivariate 
non-normal random numbers with user-specified multivariate skewness and kurtosis as 
well as the covariance matrix. If the values of the multivariate skewness and kurtosis 
are beyond the valid range of our method, the users will get a warning message and the 
allowed ranges. The package not only implements our multivariate method (function: 
mnonr), but also provides the Vale and Maurelli (1983) method (function: unonr). In 
addition, univariate and multivariate skewness and kurtosis significance tests are 
included (function: mardia). 

Example. We now illustrate how to generate non-normal data with the mnonr 
package. Suppose the goal is to generate bivariate non-normal data with multivariate 
skewness ${ = 3 and kurtosis 63 = 61. Both variables have mean 0 and variance 1. The 
covariance between them is set to be 0.5. In total, we generate 10,000 bivariate random 
numbers with the desired features. 

To generate the data, the R function mnonr was used in which we set n = 10,000, 
p= 2, ms = 3, mk = 61, and Sigma = matrix(c(1, 0.5, 0.5, 1), 2,2). The meaning of 


each argument is listed below: 


e n: the size of random number to generate; 


p: the number of variables; 


ms: the value of multivariate skewness; 


e mk: the value of multivariate kurtosis; 


Sigma: the covariance matrix of variables. 


ortoaw&ker wn 


oO 
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For illustration, we also calculated the covariance matrix of the generated data 
and conducted hypothesis testing of the univariate and multivariate skewness and 


kurtosis through the function mardia. 


The R input and output are given below. 


> mvn.data=mnonr(n=10000, p=2, ms=3, mk=61, Sigma=matrix(c(1,0.5,0.5,1) ,2,2)) 
> cov(mvn. data) 
[,1] [,2] 
[1,] 1.0795673 0.5435589 
[2,] 0.5435589 1.0378786 
> mardia(mvn.data) 
Sample size: 10000 
Number of variables: 2 


Marginal skewness and kurtosis 
Skewness SE_skew Kurtosis SE_kurt 
[1,] 0.9397589 0.02449122 24.99584 0.04897755 


[2,] 1.1900858 0.02449122 17.70874 0.04897755 


Mardia’s multivariate skewness and kurtosis 


b Z p-value 
Skewness 3.154189 5256.9822 10) 
Kurtosis 64.110259 701.3782 (0) 


The sample data yield a multivariate skewness 3.15 and multivariate kurtosis 
64.11. The covariance matrix is close to the specified one. It also shows clearly that the 
marginal univariate skewness and kurtosis for the two variables are different. According 
to marginal Equation (2), the theoretical skewness and kurtosis for marginal variables 
are: 71(%1) = 1.22, yo(41) = 29.50, y1(a@2) = 0.95, yo(a2) = 19.56. The scatter-plot and 
marginal histograms are shown in Figure 1. Even though both variables have 
leptokurtic distributions, 7, has larger kurtosis than x2, which shows on the figure that 
the distribution of x7; has a fatter tail. This is because when we form x, the 
transformation x = A€ would yield different distributions of each x;, 7 = 1,--- ,p, even 


though the €;,7 = 1,--- ,p, are tid. 


Simulation Study 


To evaluate the performance of our method, we conducted the following 
simulation study by varying the sample sizes, covariances, number of variables, and 


different combinations of multivariate skewness and kurtosis. 
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Study Design 


The sample sizes are set to be 100, 1000, and 10000. We set the variances all to be 
1 and varied the covariance between two variables from low to high with values 0, 0.1, 
0.3, 0.5, 0.7 and 0.9. In each condition, the covariances of any two variables are set to 
be the same. The numbers of variables in the multivariate data are set to be 2, 4, and 6 
which are also the number of €; in €. 

The values of multivariate skewness and kurtosis are chosen based on Cain et al. 
(2016). They provided a descriptive table of Mardia’s multivariate skewness and 
kurtosis values collected from 136 multivariate studies. We choose the minimum, first 
quartile, median, and third quartile values when the sample sizes are larger than 100. 
The values of multivariate skewness are Sf = 0, 1, 3, and 15. The multivariate kurtosis 
values are 35 = 10, 32, 61, and 91. 

We deleted some conditions due to the restricted range of multivariate skewness 
and kurtosis by Equation (16). In total, 480 conditions are evaluated, and 1,000 


replications of data are generated under each condition. 


Evaluation 


We evaluated the performance of our random number generation method by 
comparing the statistics of the generated data with the population ones used to 
generate the data. Specifically, the statistics included multivariate skewness (7, 
multivariate kurtosis 33, €;’s third moment E(£)*, and €;’s fourth moment E(€)*. We 
calculated the bias (B) and relative bias (RB) of the simulation results of the above 
statistics. The bias is the difference between the mean of the sample statistic value (0) 
and its corresponding population parameter value (@); and the relative bias is the 
proportion of the bias of the population value, which are 


ere. 


B 
N 


0, (18) 


RB = 5 x 100% (19) 
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Results 


For the sake of space, we only report several representative conditions in Table 1 
and the full results are available in the supplementary materials. They represent the 
small (67 = 1, 8} = 32), medium ((6{ = 3, 83 = 61), and large (G7 = 15, 83 = 91) 
multivariate skewness and kurtosis combinations. 

When the sample size increases, the bias of both univariate and multivariate 
measures becomes smaller. The performance of €; verified that the modified power 
method does not affect the accuracy of the power method for generating univariate 
non-normal data. For multivariate measures, kurtosis tended to be underestimated and 
skewness tended to be overestimated. Additionally, multivariate kurtosis had smaller 
relative bias than multivariate skewness. 

When comparing simulation results with various covariance settings, we found 
that both multivariate skewness and kurtosis do not seem to be affected by covariances. 
The multivariate skewness is only related to the number of variables and the value of 
E(€)°. Similarly, the multivariate kurtosis is influenced by the number of variables and 
the value of E(€)*. Covariance does not play a crucial role in multivariate skewness and 
kurtosis via our generating method, and therefore, variance-covariance matrix only 
affects the marginal measures rather other multivariate measures. 

Increasing the number of variables leads to less biased skewness and kurtosis 


estimates holding other conditions constant. 


Conclusions and Future Directions 


In this paper, we proposed a new method for generating multivariate non-normal 
data. The advantage of our method is that it allows researchers to directly specify both 
multivariate skewness and kurtosis to better control them. With the data generating 
model, we established one possible solution to relate multivariate measures to univariate 
measures: using univariate measures to generate € and apply the variance-covariance 
matrix to produce multivariate x. Our method can help researchers better understand 


the influence of the multivariate non-normality. 
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The widely used VM method and our method are both based on Fleishman’s 
polynominal related approach. Both can generate correlated multivariate random data. 
However, the two methods also have important differences. The main difference lies in 
the perspectives of the multivariate non-normality. First of all, the multivariate 
non-normality can be simply because of the non-normality of the marginal distribution 
and/or the multivariate distribution. The VM method concentrates on the univariate 
non-normality without specifically controlling the multivariate non-normality. On 
contrast, our method focuses on the multivariate non-normality but not controlling the 
univariate marginal non-normality. For instance, in a bivariate distribution of 
X = (21, 22)', through the VM method, researchers can specify the marginal univariate 
skewness and kurtosis of x; and xo, but not the multivariate measures. With our 
method, one can directly determine multivariate skewness and kurtosis of x, but with 
no control of the marginal distribution. The choice of the two methods should be based 


on the particular research interest. 


Because of the use of the power method, our method also inherits the same 
problems associated with it, such as the Gaussian-like property and multiplicity 
solution issue. First, as it is discussed by Foldnes and Grgnneberg (2015), to evaluate 
the robustness of Gaussian ML estimation using multivariate data with Gaussian-like 
property, even with the marginal univariate measures showing severe non-normality, the 
researchers might get the biased results. Without further exploration, we could not 
identify the degree of the potential impact related to our method. Second, there are 
different sets of coefficients (a, b, c, d) in the modified power method. For example, in 
the limited simulation experiment in the supplementary materials (see Table 3), we 
found that within each parameter (i.e., 6], 65, n,p) setting, there were four sets of 
possible coefficients. Different starting values will yield different sets of coefficients, 
which could affect the multivariate skewness and kurtosis. We recommend that the 
researchers should try different starting values in data generating and our R package 


provides such an option in addition to the default value. 


As shown in the simulation results, with a small sample size (n = 100), the 
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relative bias of the multivariate measures could be very high. With the increasing 
number of variables (p) and sample sizes, this issue becomes less severe. When merging 
the data of each marginal univariate variable, a small deviation could lead to a large 
gap for multivariate data. This drawback is shared with other multivariate data 
generators relating to the reliability of multivariate measures. As a future direction, we 
plan to develop a sample size planning method of different multivariate skewness and 
kurtosis to optimize the generating process. 

Since our method only used one approach to generate univariate variables, 
another related limitation as described in the method section is that some combinations 
of skewness and kurtosis cannot be obtained. However, our procedure provides 
researchers a simple data model to transform multivariate measures to univariate ones. 
In the future, we will apply other univariate generators to our method in order to 
improve the empirical performance of the multivariate generator and eliminate the 
potential problems that are related to the current modified power method such as the 


solution multiplicity and Gaussian-like property. 


Open Practices Statement: The data in this study is based on simulation. None of 


the data or materials is related to any experiments. 
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GENERATING MULTIVARIATE NON-NORMAL DATA 


Table 1 


Simulation Results (Partial) 


B(6z) B(83)  RB(Bz)% RB(B3)% B(E())  B(B(é)*) 

p=2 

(a) Bj = 1,83 = 32 

N=100 3.384 -11.147 338.4 -34.8 0.160 1.021 

N=1000 1.180 -2.853 118.0 -8.9 -0.017 -0.824 

N=10000 0.180 -0.411 18.0 -1.3 0.007 -0.159 

(b) BF = 3,83 = 61 

N=100 6.781 -29.623 226.0 -48.6 -0.032 1.238 

N=1000 3.353 -9.331 111.8 -15.3 -0.060 -2.267 

N=10000 0.654 -1.425 21.8 -2.3 0.015 -0.452 

(c)@¥ = 15, 8% = 91 

N=100 1.486 -51.942 9.9 -57.1 -0.125 -9.706 

N=1000 3.775 -17.465 25.2 -19.2 -0.081 -3.658 

N=10000 1.025 -2.816 6.8 -3.1 0.015 -0.660 
p=4 

(a)6y = 1, BS = 32 

N=100 1.066 -10.313 106.6 -32.2 -0.067 -1.377 

N=1000 0.835 -4.332 83.5 -13.5 -0.027 -0.675 

N=10000 0.227 -0.478 22.7 -1.5 0.002 -0.010 

(b) BF = 3,65 = 61 

N=100 5.834 -17.100 194.5 -28.0 -0.030 -0.761 

N=1000 1.690 -3.570 56.3 -5.9 0.018 -0.143 

N=10000 0.239 -0.239 8.0 -0.4 0.054 0.005 

(c)@¥ = 15, 83 = 91 

N=100 2.080 -36.198 13.9 -39.8 -0.062 -0.135 

N=1000 1.895 -8.388 12.6 -9.2 0.032 0.332 

N=10000 0.359 -0.805 2.4 -0.9 0.003 -0.084 
p=6 

(b) BF = 3,83 = 61 

N=100 2.554 -15.119 85.1 -24.8 -0.060 -1.265 

N=1000 1.704 -4.756 56.8 -7.8 0.023 1.077 

N=10000 0.320 -0.606 10.7 -1.0 0.011 0.048 

(c)@* = 15, 83 = 91 

N=100 2.673 -20.626 17.8 -22.7 -0.030 0.166 

N=1000 0.771 -3.770 5.1 -4.1 0.003 0.050 

N=10000 0.175 -0.157 1.2 -0.2 0.000 0.000 
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Figure 1. Scatter-plot and marginal histograms of two-variable multivariate data 
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